function [act] = find_action(X,Y,s,h)
%
% function [a] = find_action(s,h)
%
%	Guesses an optimal Q-value action given a state.
% This is the same algorithm from the updating section.

% here we try to find the optimal action for the next state by sampling
% from state-action combinations around the given state
n = 5;	% sample from a 5-degree polynomial surface
epsilon = 1e-9;	% convergence criteria
Q_next = 0;
st = s;
stp = st';
% max_iter = 10;
% iter = 1;

if size(X,1) > n
	% while iter < max_iter
		% sample n actions around the state
		% distances = sqrt(sum((repmat(st',size(X,1),1) - X(:,1:size(st,1))).^2,2));
		distances = sum((stp(ones(size(X,1),1),:) - X(:,1:size(st,1))).^2,2);
		[sorted,ix] = sort(distances);
		ix = ix(1:n);	% keep first n indices
		samples_x = X(ix,:);
		actions = samples_x(:,size(st,1)+1:end);

		% predict q-values for these actions
		for i=1:n
			samples_q(i,:) = predict(X,Y,st,actions(i,:)',h);
<<<<<<< .mine
        end
        if any(isnan(samples_q)), keyboard, end
		
		% fit a quadratic surface to the sampled points
		% A = [repmat(st',size(actions,1),1) actions];
		A = [stp(ones(size(actions,1),1),:) actions];
		A = [A A.^2];
        % tobi
        % x = A\samples_q;
        x = pinv(A)*samples_q;
	
		% find maximum of fitted function
		rhs = x(1:size(x,1)/2);
		lhs = 2*x(size(x,1)/2+1:end);
		lhs(~lhs) = ones(size(lhs(~lhs)));
		maximum = -rhs./lhs;
	
		% sample new point at maximum of fitted function
		st = maximum(1:size(st,1));
		act = maximum(size(st,1)+1:end);
		Q_pred = predict(X,Y,st,act,h);
        %if act < 0, keyboard, end
		% stopping condition	
		if abs(Q_pred - Q_next) < epsilon
			break;
=======
>>>>>>> .r89
		end

		max_i = find(max(samples_q));
		act = actions(max_i,:);
		
		% % fit a quadratic surface to the sampled points
		% % A = [repmat(st',size(actions,1),1) actions];
		% A = [stp(ones(size(actions,1),1),:) actions];
		% A = [A A.^2];
		% x = A\samples_q;
		% 	
		% % find maximum of fitted function
		% rhs = x(1:size(x,1)/2);
		% lhs = 2*x(size(x,1)/2+1:end);
		% lhs(~lhs) = ones(size(lhs(~lhs))); 
		% maximum = -rhs./lhs;
		% 	
		% % sample new point at maximum of fitted function
		% st = maximum(1:size(st,1));
		% act = maximum(size(st,1)+1:end);
		% Q_pred = predict(X,Y,st,act,h);
		% 	
		% % stopping condition	
		% if abs(Q_pred - Q_next) < epsilon
		% 	break;
		% end
	% 	Q_next = Q_pred;
	% 	iter = iter + 1;
	% end
else
	% don't know value
	act = [];
end
